Generation of dispersive shock waves by the flow of a Bose-Einstein condensate past a 

narrow obstacle 
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We study the flow of a one-dimensional Bose-Einstein condensate incident onto a narrow obstacle. 
We consider a configuration in which a dispersive shock is formed and propagates upstream away 
from the obstacle while the downstream flow reaches a supersonic velocity, generating a sonic horizon. 
Conditions for obtaining this regime are explicitly derived and the accuracy of our analytical results 
is conflrmed by numerical simulations. 
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I. INTRODUCTION 



The flow of a quantum fluid past an obstacle can gene- 
rate various excitations. In the classical example of liquid 
Hell these arc phonons and rotons introduced by Landau 
[l| in his theory of superfluidity. According to Landau, 
superfluidity is lost when the flow velocity past an obsta- 
cle or through a capillary tube exceeds the threshold of 
creation of Cherenkov-like radiation of linear waves (or, 
in other words, of quantum quasi-particles). In actual 
experiments where Hell flows through a narrow chan- 
nel, superfluidity is lost at velocities much lower than 
what is predicted by the Landau criterium (see, e.g., 0). 
This discrepancy between theory and experiment was ex- 
plained by Feynman Q as resulting from the nucleation 
of nonlinear excitations: vortex rings generated by the 
flow past the obstacle. 

The same phenomenology is observed in the flow of 
Bose-Einstein condensates (BEC), either for trapped ul- 
tracold atomic vapors or polaritons in semiconductor mi- 
crocavities, with special features linked to the specific 
dispersion relation of elementary excitations in these sys- 
tems: for instance the Landau critical velocity is the 
velocity of sound c, and in two dimensions (2D) the 
Cherenkov radiation forms an interference pattern lo- 
cated outside of the Mach cone 043 ■ These specifici- 
ties being taken into account, one observes phenomena 
similar to those observed in liquid helium: in 2D or 3D 
flows past obstacles, superfluidity is broken at velocity 
lower than c [1, because of the nucleation of vortices 



IC - 14 1 or generation of effectively stable oblique solitons 
15l-ll7l| (recently observed in experiments with polariton 
condensates [H, E^), and more complicated dispersive 
shock waves (DSW) patterns ^2^^. 

The situation in quasi- ID flows is similar: although 
generation of vortices or oblique solitons is impossible 
in this case, dark solitons are still easily generated stable 
nonlinear excitations 2^ ■ Together with DSW, these 
nonlinear excitations are, as is already the case in higher 
dimension, keys ingredients for understanding the time- 
dependent dynamics of guided BECs in presence of obsta- 
cles, as experimentally studied in Refs. p5l - [27| . Besides 
their interest for studying nonlinear quantum transport, 
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FIG. 1. (Color online) Density profile in vicinity of the delta 
potential represented by a (red) vertical straight line. The 
flow is stationary, with a velocity directed toward positive x 
as indicated by the arrow. The density in the region i < is 
a portion of a dark soliton (see Sec. IIII B|l with asymptotic 
upstream density (ni) and subsonic velocity. The flow in 
the downstream region a; > has a constant density and a 
supersonic velocity. 



quasi-lD BEC flows have also been suggested as model 
systems for creating sonic horizons suitable for the exper- 
imental observation of acoustic analogues of Hawking ra- 
diation [28l - [33 |. In particular, an interesting type of sta- 
tionary sonic horizon has been identified in Refs. [33.[35|. 
It corresponds to the situation where an upstream sub- 
sonic region is separated from a downstream supersonic 
one by an obstacle of the form of a delta potential (see 
Fig-ffland discussion below). The present work is devoted 
to the study of the dynamics of the formation of this con- 
figuration. We show that it can be obtained by launching 
a BEC wave packet onto a localized obstacle (not neces- 
sarily a delta peak [36']). When the wave packet reaches 
the potential, the density typically piles up in front of 
the obstacle, forming a plateau accompanied by a DSW 
which is ejected upstream. We study the characteris- 
tic features of this DSW in detail (both analytically and 
numerically) and determine for which parameters (spe- 
cific to the flow and to the obstacle) the situation just 
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described can be realized. 

Theoretical analysis of ID flows past an obstacle has 
already been addressed in a number of papers (see, e.g., 
[35., .37-40J). In these papers, it was found that there 
exist two critical velocities U- and u+ (u- < c < 
The sub-critical flows whose velocity u is below u_ are 
superfluid and generate no excitations in vicinity of the 
obstacle; the super-critical ones with velocity above Uj^. 
do not generate nonlinear DSW but the Cherenkov ra- 
diation is effective and there is no sonic horizon; at last, 
in the trans-critical region, for flows whose velocity u lies 
between the two critical values, < u < u+, DSW are 
generated, generally speaking at both sides of the obsta- 
cle, and only in a very narrow range of velocities both 
DSW are detached from the obstacle. Hence, it seems 
difficult to reach a situation such as the one illustrated 
in Fig. m 

Note however that the above quoted references [35l.[37l- 
|40| mainly considered flows with identical asymptotic pa- 
rameters at both sides far enough from the obstacle. This 
choice of boundary conditions is natural when the obsta- 
cle is put into motion starting from a situation where it 
is immobile in a condensate at rest. However, another 
setup is possible, which is of considerable practical inter- 
est, being similar to the flow of a fluid through a Laval 
nozzle [4ll|. In this configuration (already considered in 
Ref . (44I ) the parameters of the flow are kept fixed only 
on the side of the incoming flow and the downstream flow 
expands freely into vacuum. In this case a DSW can be 
created only upstream and the parameters of the shock 
can be calculated analytically for two typical models of 
obstacle potential: (i) a smooth potential with typical 
size greater than the characteristic healing length ^ and 
(ii) a short-range potential acting at distances much less 
than ^. In the flrst case the dispersion effects can be ne- 
glected at the location of the obstacle (in the so called 
hydraulic approximation) and the theory reduces to the 
scheme presented in Ref. [i^. In the second case the 
potential can be approximated by a J function and its 
action amounts to the matching condition of exact solu- 
tions at the point of its location. This last case has not 
been considered so far in this kind of problem and we 
shall discuss it here in detail. We will show that it makes 
it possible to realize asymptotically (i.e., at large time) a 
sonic horizon such as represented in Fig. [TJ 

The paper is organized as follows. In Sec. |II]we present 
the problem and the typical dynamical situation we aim 
at describing. In Sec. IIIII we discuss the time-dependent 
analytical solutions of the flow in each of the charac- 
teristics regions of space identified in the previous sec- 
tion. In Sec. IIVI we briefly compare our results with the 
ones obtained in the case of a thick obstacle. In Sec. 
IVl we compare the analytical results with numerical si- 
mulations and propose an improvement of the analytical 
description. Finally we present our conclusions in Sec. 

lYIl 



II. FORMULATION OF THE PROBLEM 

In the so-called ID-mean field regime [1^ the dynamics 
of the condensate is described by the Gross-Pitaevskii 
equation which governs the evolution of the wave func- 
tion ^{x,t). Expressing densities in units of a refe- 
rence density nref, energies, distances and velocities in 
units of the corresponding chemical potential /ircf (f^rcf ), 
healing length ^rcf = fi-/(™Mrcf)^^^ and speed of sound 
Crof = fi / (fn^rcf) , one can write the Gross-Pitaevskii 
equation in the form 

ii^t^-^i'.cc + [uix) + \^p\'']ip . (1) 

The Gross-Pitaevskii equation is a sufficient ingredient 
for describing the generation of DSW in a BEG as shown 
by the comparisons between theory and experiments dis- 
played in Ref s. [13, m. 

By means of the Madelung substitution 



ij{x,t) = y^n{x, t) exp l^i j u{x',t)dx'j e"'^* (2) 

the Gross-Pitaevskii equation can be cast into an 
hydrodynamic-like form for the density n{x,t) and the 
fiow velocity u{x,t): 

rit + {nu)x = , 

/ n1 nxx\ rr (3) 

ut + uux + nx+ \ — = -Ux ■ 

We now briefly introduce the concept of Riemann in- 
variant by considering the very simple example of the 
dispersionless limit of Eqs. ([3]). This case is not treated 
only for pedagogical purposes. At the boundary between 
the DSW and regions of flat profiles (where dispersion is 
indeed negligible) it will make it possible to match the 
description of the non linear wave in terms of elaborate 
Riemann variables [Eqs. (jTS]) and p4)) ] with a simple 
description of the type specified by Eqs. ([5]), ^ and ([7]) 
[this matching will be achieved by means of Eqs. ([25)1 
and (P5|)]. 

The last term of the left hand-side of the last of Eqs. ([3]) 
is responsible for the dispersive character of the BEG 
wave. In the absence of potential, and in a regime where 
the effects of dispersion can be neglected, Eqs. ^ reduce 
to 

nt + {nu)x = , tit -f UU.J, -f n^, = . (4) 

These equations can be written in a more symmetric 
form by introducing the following Riemann invariants 



A±(x,t) = ^^±y;^, (5) 
which evolve according to the equations [equivalent to 

9tA±+i;±(A+,A_)5,A± =0, (6) 
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with 



^^±(A+,A_) = i(3A± + A, 



(7) 



We will encounter below other Riemann invariants which 
describe the DSW and obey equations similar to ^ and 
([7]) [Eqs. (ITOl) and d^ ]. However, the dispersive nature 
of the shock will be there fully taken into account, in 
contrast to the simple approximation ^ and ([7]). 

Let us now present the initial velocity and density pro- 
file of the incident flow. We suppose that at initial time 
(t = 0) a half-infinite pulse of BEC with a step-like dis- 
tribution 



iP{x,t = 0) 







no expliwox} for x < , , , 
for x>0 . 



collides with a repulsive 5 potential located at the origin 
of the coordinate system; 

U{x) = Kd{x), K>0. (9) 

In other words, the initial density and flow velocity are 

riQ for a; < , 



n{x, t — 0) = 



for a- > , 



u(x, t = 0) = 



uo{> 0) 




for 
for 



X < 
x > 



(10) 



(11) 



At later times, the typical situation we aim at describing 
is illustrated in Fig. [2] which represents the density dis- 
tribution n{x, t) at some fixed time t > 0. It corresponds 
to a flow initiated by the profile pH)) and (ITT]) in which 
a plateau develops upstream of the potential (region C) 
while a dispersive shock wave (region B) propagates away 
from the obstacle, in the negative direction. The flow just 
downstream from the obstacle forms a supersonic plateau 
(region D). The right edge of region D matches with a 
simple- wave solution (not shown in Fig. which descri- 
bes how the density vanishes at large x. In the present 
work we are not interested in the description of this part 
of the flow pattern. As one can see, the typical fiow we 
consider can be subdivided in this case into several re- 
gions (denoted as A, B, C and D in Fig. ^ with specific 
features in each region. 

• In region A [x < X-{t)] we have the incoming flow 
ipA with the parameters (|10l) and PT|) . This flow can be 
considered as a boundary condition: 



nA{x,t) = no , 
UA{x,t) = Uq , 



X < X^(t) , for aU t > 0. (12) 



In this region there is no dispersion nor external potential 
and the analysis in terms of Riemann invariants ([5]) is 
trivially valid: the Riemann invariants are constant with 

• In region B [X^{t) < x < X^{t)] a dispersive shock 
wave takes place which can be described as a modulated 



n(x,t) 




XJt) XJt) 



FIG. 2. (Color online) Sketch of the density profile n{x,t) at 
a given time t > 0. The flow is directed toward positive x. 
The regions denoted as A, B, C and D are presented in the 
text. Points X-{t) and X+{t) are the small amplitude edge 
and soliton edge of the DSW (region B). The characteristic 
densities no, ni and n2 are defined in the text [Eqs. (llOp . (|24p 
and (l26l)]. 



periodic solution ipnix, t) of Eqs. ([T]) and ^ with U(x) = 
0. Such a solution can be written in the form (see, e.g., 
0,0) 

nn{x,t) =i(A4-A3- A2 + Ai)2+ 
(A4 - A3)(A2 - Ai)x 

sn2 (V(A4 - A2)(A3 - Ai) {x - Vt),m}j , 

(13) 



UB{x,t) = V - 



C 



(14) 



nB{x,t) 

where sn is the sine elliptic Jacobi function, 

(A4- A2)(A3-Ai) 

and 

C = i(-Ai - A2 + A3 + A4) X (-Ai + A2 - A3 + A4)x 



(Ai - A2 - A3 + A4) 



(16) 



In strictly periodic solutions the parameters Ai < A2 < 
A2 < A4 are constant and they determine characteristics 
of the wave such as the amplitude of oscillations 



a = (A2 - Ai)(A4 - A3) 
and the wavelength 

2K(to) 



L = 



V(A4- A2)(A3-Ai) 



(17) 



(18) 



K(m) being the complete elliptic integral of the first kind. 
In the limit m — )■ (A2 = Ai), describes a small 
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amplitude sinusoidal wave, and in the opposite case m — 
1 (A2 = A3) it describes a dark soliton. 

In the case of a slowly modulated dispersive shock wave 
such as occurring in region B, the A's are functions of 
X and t which vary weakly over one wavelength and one 
period. Their slow evolution is governed by the Whitham 
equations (ieHiTj 



OK 
dt 



+ Wi(Ai,A2,A3,A4) 



dx 



0, i = 1,2,3,4. (19) 



Comparing with Eqs. ([6|) and ([7]) one sees that the A^'s 
are the Riemann invariants of the Whitham equations. 
The ViS are the characteristic velocities. Their A^'s de- 
pendence is much more complicate than the simple linear 
combinations appearing in ([7]); they can be expressed in 
terms of complete elliptic integrals of the first and the 
second kind [H, [4^ . One can use the following conve- 
nient formula for their computation [i^, H^] 



2dL/dh' 



1,2,3,4 



(20) 



At X = X- (t) we have the "small amplitude edge" of the 
dispersive shock wave [with X2{X-{t),t) = Xi{X-{t),t)] 
where the wave should satisfy the matching conditions 
with the flow in the region A. This implies that the 
mean values of the density and the flow velocity 
coincide with uq and mq respectively: 



ns{X^(t),t) = no, (X_ (t ) , i ) = wq . 



(21) 



Since at the left edge of the DSW we have A2 = Ai the 
Whitham equations for A3 and A4 can be shown to 
simplify to 



dtX3+^ (3A3 + A4) d^Xs - 
dtXi+^ (A3 + 3A4) a,A4 = 



(22) 



and these equations can be identified with the Riemann 
form ^ and ([7]) of the dispersionless limit ^ of Eqs. (U) 
(without potential); i.e., one has 



Xi{X^{t),t) = Xl = ^uo + ^=Xo 
X3{X^{t),t) = X± = ^uo- . 



(23) 



The other edge of the DSW occurs a,t x = X+{t). We 
denote it as the "soliton edge" because at this point 
the density oscillations are soliton-like : X2{X+{t),t) — 
A3(X_|_(t), t) (i.e., m = 1). The matching conditions at 
X+ read 



nB{X+{t),t) = ni , UB{X+{t),t) = ui 



(24) 



where we suppose that X+ is located far enough from the 
origin in order that the stationary solution ipc in region 
B reaches its asymptotic values ndx) — > ni, Uc{x) — >■ ui 
when X ~ X^ (i-e., in the formal limit x ^ — 00 in region 
B, see below for more details). 



Since at the soliton edge A2 = A3, one can show as 
previously that the Whitham equations for A4 and Ai 
reduce to a form similar to ^ and ([7]) where A4 plays 
the role of A+ and Ai plays the role of A_. Hence one 
has 



Xi{X+{t),t) = A^j: = iui + V"i 
Xi{X+{t),t) = X'i = \ui-^i. 



(25) 



• The region C corresponds to a smooth and statio- 
nary solution ipc{x,t) of the Gross-Pitaevskii equation 
for X+ < X < and the region Z) to a smooth and 
stationary solution ipo for < x < Xgw The flow at 
X > Xsw represents a simple wave solution (with small 
dispersive corrections) connecting region D with vacuum 
and is of no interest to us (this is the reason why we do 
not show the region x > Xgw in Fig. [2]) . In region D we 
have a uniform flow 



no{x,t) = n2 = const 
Uo{x,t) ^ U2 = const 



for < X < X^ 



(26) 



Since the flow is stationary in both regions C and D, the 
parameters (|24p and (|26p must satisfy the condition of 
conservation of flux 



(27) 



and the Bernoulli law (constant value of the chemical 
potential /i) 



1„,2 



+ ni = 2^2 + n2 = n 



(28) 



One must also satisfy the condition of continuity of the 
wave function at x = 0, 



(29) 



and the jump condition for the derivative of the wave 
function which follows from Eqs. ([1]) and © 



9,1^0 (0, i) - 9,i^c (0, i) = 2 K (0, i). 



(30) 



Thus, our task is to find the solution which satisfies all 
the above matching conditions [Eqs. (gTl), ([Ml), (HZl), (E3), 
and (|5n|) ] and yields the parameters rti, mi, ri2, U2 
and the positions X± {t) of the edges of the DSW as func- 
tions of the incoming fiow parameters rip, uq and of the 
potential strength k. 

It is clear that such a solution does not exist for any 
choice of the parameters no, uq, k. For example, in the 
limit K — >■ our system reduces to the well-known "dam 
problem" 0,|43| which is described in the hydrodynamic 
approximation by a simple wave solution without forma- 
tion of dispersive shock wave. Therefore our solution 
should be complemented by an explicit statement of its 
conditions of existence. 
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III. SOLUTION OF THE PROBLEM 

A. Dispersive shock wave (region B) 

As was indicated above, the DSW is described by four 
parameters A^, i = 1, 2, 3, 4, which change slowly along 
the wave. We suppose that the DSW is detached from the 
obstacle and propagates upstream with the velocities V± 
of the edge points X± . The distance at which tpc reaches 
its asymptotic value (with ric — ni, Uc — ui) when one 
goes away from the origin is of order of the healing length 
in the region C, i.e. ~ n-^ Since this distance is 
much less than the (time increasing) length \X^\, we can 
assume in good approximation that the formation of the 
plateau ric — ni, Uc — ui occurs instantaneously. Then 
the DSW can be described by a self-similar solution of 
the Whitham equations (|19p and for the left-propagating 
DSW we get 



Ai = const 



A., 



const 



A. 



const 



and 



W2(Ai,A2(x/i),A3, A4) = - 



(31) 



(32) 



1 - 



-1 - 



t 






region 


region B 


region 


A 




C 



















x/t 



FIG. 3. (Color online) The Riemann invariants plotted as 
functions of the self-similar variable x/t in regions A, B 
and C. In region C, the (red) dashed line (not consid- 
ered in Section IIII|) corresponds to the soliton train with 
A2 = A3 = x/t — 111/2 discussed in Section |V] [Eq. (f86|) ]. 
The figure is drawn for the situation studied in Sec. [V] with 
no = 1, Mo = 1 and n — 5.2. The corresponding veloci- 
ties of the edges of the DSW (region B) are V- = —2.4 and 
V+ = -0.48. 



This is the so-called Gurevich-Pitaevskii problem intro- 
duced into the theory of DSW in Ref. [51| and used for 
description of internal waves generated by the flow of wa- 
ter past an uneven bottom in [H, [111 (see also [i^l and 
references therein). 

When ni and ui are found, the values of the constants 
Ai, A3 and A4 can be determined from Eqs. (^5)) and (^5)) 
whereas Eqs. (^0]) and determine A2 as a function of 
x/t: 



V2(Al, A2, A3, A4) 



A,+ 



(A3 - A2)(A2 - Ai)K(m) ^ X 
(A3 - A2)K(m) - (A3 - Ai)E(m) t 



(33) 



The plots of Ai, i — 1,2,3,4, as functions of x/t are 
displayed in Fig. [31 The left edge of the DSW moves 
with velocity 



V- 



= w2(Ai, Ai, A3, A4) 



= i(2Ai + A3 + A4) + 



2(A3- Ai)(A4-Ai) 



(34) 



2Ai — A3 — A4 
and the right edge moves with velocity 

V+ = ^= V2{Xi,X3, A3, A4) = i(Ai + 2A3 + A4). (35) 

For finding ni, ui and, hence, Ai, we have to turn to the 
solution of the Gross-Pitaevskii equation in region C and 
its matching with the solution in region D. 



B. Flow across the S potential (regions C and D) 

In regions C and D the flow is stationary, tp{x,t) — 
e~'^*(^(a;), so that the Gross-Pitaevskii equation reduces 
to 

- l^xx + |<^P<y5 -I- K6{x)ip = (iuj -I- ni) ip. (36) 
It is convenient to introduce temporarily the variables 



(37) 



Mc being the asymptotic (formally at x — > —00) Mach 
number of the flow in region C. Then Eq. (|36l) takes the 
form 

- \cj)yy + + KSiy)^ = (iM2 + 1) 0. (38) 
From conservation of the flow (|27|) we get 

U2 _ Tlx 

ui n2 

and 77 can be found as a function of Mc — u\/ ^/n\ < 1 
with the use of Eq. ([28l) or 



(39) 
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1 + ^) '7 + 1 = 0. 



(40) 



The relevant solution (77 > 1) of this equation is given by 
1 ' ' 



7] 



(41) 
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We look for a solution of Eq. (|38| of the form 



C. The global solution 



0(2/) 



— JyMc 



{cos 9 tanh[(y — ijq) cos 9] — i sin 9} , (42) 



for y < and 



(43) 



for y > 0. In the above expressions (|42l) and P5|) we have 
unl \fn\ = ?7Mc and sin0 = Afc- This solution has been 
first identified in Ref. [sl]. The flow upstream from the 
obstacle corresponds to a portion of a dark soliton which 
is attached at j/ = to a downstream supersonic fiow. 
The condition of continuity of the function 0(?/) at y = 
[see ([29])] gives 

sin 7 = y/r] sin 9, cos 7 = tanh(yo cos 6*) , (44) 

from which we get 



tanh(yo cos ( 



1 - ryM2 



(45) 



V '/(I - ) ■ 

The condition ([50)) takes the form 

</'!;(O+)-0y(O") = 2«0(O). (46) 



Substitution of expressions (|42|) and (143)) in this relation 
gives after some algebra with the use of Eqs. ()44l) and 
(gSj) the equation 



(47) 



Elimination of 77 with the help of Eqs. (|4T)) yields 

?i = F(M,) i.e., K = V^^^(Me), (48) 

where 

3/2 

F(M,) = :^ I Jl + --,-3l . (49) 




Hence, according to (|57)) 



F(Me 



(50) 



Then Eqs. (HHj) and dH]) yield the expressions for n-i and 
U2: 



W2 = Ml?? 




(51) 



(52) 



We thus have obtained the parameters ni, mi, 712, ui as 
functions of k and M^- We now have to relate to 
the physical parameters rig and uq which describe the 
incoming flow. 



To get the global solution, we use the relation = 
A4 = \% [see Eqs. ^ and ^ and Fig. [5] and (gS]) and 
(iroi) to obtain 



Ao = + 2"i + = YijA^ \ 

(53) 

As a result, all expressions can be written in a parametric 
form with playing the role of the parameter. From 
(|53)) we get 



k(Ao, Afc) = Ao 



1 + Afc/2 



(54) 



which determines Mc as a function of k and Aq. Substi- 
tution of this expression into (|48 )) - (|50)) yields 



(1 + A'4/2)2 



1 + Afc/2 ' 
(55) 



n2(Ao,Afc 



U2(Ao, A#c) = 




(57) 



In combination with (|54)) . formulae (|55)) . (|56)) and (|57)) 
determine - in a parametric form - the dependence of 
ni, wi, 712, U2 on K and Aq = ^uo + v^^^o- These formulae 
represent an important step in our study since they de- 
termine the overall structure of the flow (ni, ui, 712 and 
M2) as a function of the initial parameters characterizing 
the incident beam (779 and mq) and the obstacle (k) with 
which it collides. 

Actually, the overall structure of the fiow is fully cha- 
racterized only once the locations X^(t) and Xj^(t) of 
the boundaries between the different zones are known 
(cf. Fig. [2]). According to our self-similarity hypothesis 
X±{t) = V± t, where V± are the time-independent veloci- 
ties of the boundaries. Now that 77i and ui are found, V+ 
and V- are simply determined as follows. One first can 
get [from Eqs. ()23p and (|25)) ] the expressions for the 3 
Ricmann invariants which remain constant in the region 
of the DSW 



Ai 



-Ao 



2~M, 
2 + Mc 
A4 



A3 



(58) 



Then, substitution of formulae (|55)) into Eqs. (|M)) and 
(|35)) yields the desired expressions of the velocities V- 
and V+ of the edge points of the shock wave. In par- 
ticular, it is interesting to realize [cf. Eq. dSTt ] that the 
velocity V- of the small amplitude edge of the DSW 
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corresponds to the group velocity of linear excitations in 
region A of the system with wavelength L{X_), where L 
is the (position dependent) wavelength of the nonlinear 
oscillations ([T^. This is obtained by noticing that the 
explicit formula (p4)) and the expressions ((58)) of the Rie- 
mann invariants makes it possible to write V- under the 
form 



1 + Mc/2 



.(2- 



- 2 - 2Mo 
Mo){Mo-M,) 



(1 + M,/2y 



(59) 



where Mq = uo/y^no. Similarly, when Ai — A2, the ex- 
pression (|18l) for L reduces to 



L{X_ 



V(A4-Ai)(A3-Ai) 
7: {1 + MJ2) 

y^{Mo + 2){Mo- Mc) 



(60) 



Defining k = 271 ^/fiQ / L{X ^) as the dimensionless wave 
vector at the small amplitude edge of the DSW, it is 
a simple exercise to show that using Eq. (|60p one can 
rewrite (|59l) under the form 



= Mo - Co 



1 + fcV2 



(61) 



which is the group velocity of excitations whose disper- 
sion relation is given by a;(fc) = u^k — c^k^Jl + 
where cq — i/uq is the speed of sound in a uniform sys- 
tem at rest with density uq. This dispersion relation 
corresponds to (left propagating) linear elementary exci- 
tations in a BEC with constant density uq and constant 
velocity uq. Relation (pT|) is very natural since V_ is the 
velocity of the small amplitude edge X- at which the 
DSW is linear: the linear front of a wave-packet should 
propagate with the appropriate group velocity. This is 
a generic feature of the small amplitude edge of a DSW, 
as proven on general grounds in Ref. [s^l- Relation ([6T|) 
nonetheless constitutes a non-trivial check of the validity 
of our description, since the connection between V- and 
L{X^) is not straightforward. 

As for the soliton edge of the DSW, we obtain from 
Eqs. §SI and (E 



(62) 



V+ 



tto(l + M,)-27?i^ 
2 + Mc 



This is the velocity of the edge soliton whose amplitude 
is equal to 

a = (A3 - Ai)(A4 - A3) = -^^ ^^^^ . (63) 

Note that our approach not only determines the gross 
features of the flow (those which are depicted in Fig. [21 
Ml, rii, U2, n2 and X±) but also yields a precise predic- 
tion for the density profile in the region of the obstacle 



and of the DSW. In the region of the obstacle the order 
parameter is determined by Eqs. and In the 

region of the DSW, the velocity and the density profile 
are determined by Eqs. and ([H]) which need the 

input of the A^'s [Eqs. ([32]) and (|58)) ]. These predictions 
will be tested against numerical simulations in Sec. [V] 

The solution found here is based on two assumptions: 
(i) that the discontinuity ni > no arises after collision of 
the BEC pulse with the obstacle. This is equivalent to 
the condition a > or 



Uq 



(64) 



and (h) that the DSW is detached from the obstacle, i.e., 
V+ <0or 



Up ^ 2 

\/rio 1 



Mr. 



(65) 



If the density uq is fixed, then Eq. ([M]) gives the lower 
velocity uq below which - and Eq. (|65p gives the upper 
velocity above which - the solution we are interested in 
disappears. These critical values of uq are functions of 
the incoming density uq and of the potential strength k 
determined by the equations 



Uq 



and 

Uq 



F 



F(2X!^-1 

Mo 



(lower boundary), (66) 



(upper boundary), (67) 



where we recall that function F is defined by . These 
two curves are plotted in Fig. S) 

For small nj ^Juq we get the series expansions: 



_^ 3 / K 



2 \ Jn^ 



2/3 



, (lower boundary), (68) 



and 



Mo 



3 / K 

1 + T 



2/3 



, (upper boundary). (69) 



/no 4 V\/^. 

In the opposite limit of large k,/ y/np we get 

Mo \ I K, 



^/rio 2V2 Ky/no 



, (lower boundary), (70) 



and 



"° -i 2 - 2^2 ( ] , (upper boundary). (71) 
/no V Jn()J 



We emphasize here that Eqs. (|66)) and (|67)) are im- 
portant results of our study because they determine the 
region of parameters for which our analytical solution ex- 
ists, i.e., for which (i) the DSW is expelled upstream from 
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1.5 



1 - 



0.5 



K Nn, 



FIG. 4. Boundaries of the region of parameters for which the 
DSW is expelled upstream from the obstacle. The vertical 
axis is the speed of the incident beam in units of the speed 
of sound Co = y/no- The horizontal axis is the dimensionless 
strength K/^/no of the 5 potential. The DSW is expelled from 
the region of the obstacle when the point representative of the 
system lies between the two solid curves which correspond to 
Eqs. (|66|) and (|67|) . The dashed curves are the approximate 
boundaries ^E^, (ED) (valid when k/^^ < 1), and (ffT]) 
(valid when n/^/no 2> 1). 



the obstacle and (ii) the flow forms in vicinity of the ob- 
stacle a sonic horizon such as schematically represented 
in Fig. [H 

In this line, it is important to notice that the above ap- 
proach is not applicable for positive values of the soliton 
edge velocity V+ > 0, i.e., when the soliton edge of the 
DSW is attached to the obstacle. This is different from 
what happens in the case of thick obstacles [40] where 
the existence of a characteristic length I of the poten- 
tial representing the obstacle plays a crucial role. Since 
in this case I is large compared with the wavelength of 
the (possibly attached) DSW, the regions of DSW and 
of the "hydraulic solution" in vicinity of the obstacle are 
well separated and one can safely assume formation of 
a plateau at the left boundary of the region of the hy- 
draulic solution. In the present case of a (5 potential such 
a plateau forms at a distance of the same order of ma- 
gnitude as the DSW wavelength and, hence, in the case 
where the DSW is attached to the obstacle the solution 
of region C (see Fig. [5]) completely disappears and the 
formulae derived above loose their meaning. 



IV. COMPARISON WITH THE FLOW PAST A 
WIDE PENETRABLE BARRIER 



In this section we briefly compare the results we ob- 
tained for a thin potential with those obtained in Ref. 
poj for a wide and smooth potential U{x) which takes 



its maximal value at a; = 0, 

[/,„ = max{f/(x)} = U{0), (72) 
and differs from zero only inside the region 

-l<x<l- (73) 

If I and all lengths characteristic of the potential are 
much larger than the healing length of the conden- 
sate, the transition from an upstream subsonic flow to a 
downstream supersonic one is described by the so-called 
"trans-critical flow" [Z^]. It is obtained for an incoming 
velocity uq G where the critical velocities u± 

are the roots of the equation 



! = [/„ 



(74) 



(to simplify the notation, in this Section we assume 
n-ici = '^o)- If f^m li the critical velocities are given by 
the series expansions 



u± 



1± 



3f/„ 



(75) 



Flow velocities ui,2 at the boundaries of the hydraulic 
region are roots of the equation 



u 

T 



Uq 



Uq 



2/3 



(76) 

The smaller root corresponds to the upstream velocity 
Ml and the larger one to the downstream velocity U2. If 
Um < 1 we get 



ui = l + i(ito - 1) - 

U2 = l+ ^{Uo - 1) + 



2U„ 



2U„ 



(77) 



The DSWs can be attached to the trans-critical flow as 
accounted for in detail in [40| . In our initial value prob- 
lem, with a wave packet incoming from the left infinity, 
there can be no downstream DSW. However, an upstream 
DSW exists. It is detached from the obstacle only if the 
condition uq < 2 — ui is fulfilled. Thus, the solution we 
are interested in, with a DSW expelled upstream from 
the region of the obstacle, corresponds to a region in the 
plane (J7,n,ito) where 



U_ < Uq < 2 — Ui . 



(78) 



This determines the parameters of the flow and of the 
potential leading to the formation of an acoustic analog 
of a black hole in the flow of a condensate past a thick 
obstacle. If [/m ^ 1, then inequalities (l78t take the form 



3Un 



<Uq <1 + 



3Un 



(79) 



The whole region (1751) is shown in Fig. [5] which is the 
thick obstacle analog of Fig. 
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FIG. 5. Boundaries of the region (|78|) for which the hydrauhc 
solution exists and the upstream DSW is detached from a 
thick and smooth obstacle (here we use units such that no — 
1). The dashed lines correspond to the approximations (|79p 
which are valid when Um <C 1. 



V. COMPARISON WITH NUMERICAL 
SIMULATIONS 

Our analytical theory is based on the assumption that 
the stationary flow with a plateau in region C is formed 
instantaneously after collision of the condensate with the 
obstacle. This assumption justifies the application of the 
self-similar solution [characterized by formulae (|3T|) and 
([5^ ] which only depends on the variable x/t. However, 
in practice (i) any initial wave-packet has — contrarily to 
our simplified initial flow ^ — a finite region of transi- 
tion from a constant plateau density to vacuum, and (ii) 
the stationary solution in region C forms in a period of 
time negligibly small compared with the asymptotic val- 
ues of time under consideration, but nonetheless finite, 
and this may induce some noticeable effects. Therefore 
we performed numerical simulations for determining the 
limitations of our analytical approach. 

In our numerics the incident condensate wave packet 
was represented by the following initial density distribu- 
tion: 



n{x, t = 0) = no 



1 - tanh(x/A) 



(80) 



According to the prescription of Sec. |TI] we take here 
JT-rcf = ?iO: so that riQ = 1 henceforth. The velocity of the 
incident flow it taken to be uq = 1 (i.e., uq = cq = Crcf). 
We consider three different cases with different thick- 
nesses A of the incident front of the wave packet: A = 20, 
10 and 5 (in units of — 6of)- The obstacle is modeled 
by a Gaussian potential 



U{x) = Uo exp{-a;Vcr^} 



(81) 



The results of the numerical evolution of the flow after 
a time t — 500 are shown in Fig. [SI where the different 
numerical curves correspond to different values of A. The 
potential (l5Tt deviates significantly from a S potential: 
we have here a — 0.5 which is too large for using the 
natural prescription for defining, starting from (|8ip . an 
equivalent S potential of the form 1^ hy k — J U {x)dx = 
Uoa^/n = 3.5. This prescription would be accurate only 
in the limit cr <C 1. The value k = 5.2 used in the 
analytical procedure results from a fit of the numerically 
obtained density of the plateau {ni — 2.192); then all 
the other parameters of the flow are determined by the 
analytical formulae. That this procedure is legitimate is 
confirmed by the very good description of the upstream 
velocity ui in region C and of the density n2 and velocity 
U2 in the supersonic downstream region D: one obtains 
analytically ui = 0.0390, ns = 0.0412 and U2 = 2.075, 
whereas one finds numerically ui = 0.0388, n2 = 0.0410 
and U2 = 2.075 (see also Fig. [7]). 

We checked that the natural prescription k = V^a^pK 
yields a very good agreement of the analytical and of the 
numerical approaches for low values of a. For instance 
when C/q = 4 and cr = 0.1 the natural prescription yields 
K = 0.71 and for this value of k one obtains analytically 
Tlx — 1.658, whereas one find numerically tii — 1.6703. 
In this case perfect agreement with the numerical va- 
lues of ni, 712, u\ and U2 is found by using an effective 
K — 0.7365, close to the natural prescription. However, 
for reducing the numerical effort we performed extensive 
simulations only in the case a = 0.5. The efficiency of 
the effective 5 barrier for representing a potential of the 
form (|8ip even though a is not very small is a positive 
test of the robustness of our analytical description of the 
system. 

We note that the theory predicts the conservation of 
the Riemann invariant 



A^(= Ao = 1.5) = = Ml/2 + 



(82) 



with /7o = 4 and a = 0.5 (i.e., the size of the obstacle is 
less than the healing length). 



across the DSW (cf. Fig. O and that this is fulfilled 
with very good accuracy in the numerical simulations. 
The positions X±{t) of the edges of the DSW are also in 
good agreement with the theoretical prediction. Because 
of the scale used on the a;-axis of Fig. [6l one could think 
that the density is discontinuous around a; = 0. This 
is not the case, and furthermore one can verify that in 
this region the analytical and numerical density profiles 
are very similar, as shown in Fig. [7] (right plot). Also 
inside the DSW, the nonlinear oscillations are very well 
described by the analytical approach (same figure, left 
plot). 

Hence, we can legitimately state that the hypothesis of 
rapid formation of the DSW which, as already explained, 
is at the heart of our assumption of self-similarity, is va- 
lidated by the comparison with numerical simulations. 
This makes it possible to bypass the delicate and inte- 
resting question of the short time dynamics which is dis- 
cussed for instance in Ref. [4^ . However, our numerical 
simulations reveal a peculiar feature of the flow, namely. 
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FIG. 7. (Color online) Density profiles at t = 500. These 
plots are enlargements of specific parts of the two lower plots 
of Fig. H The left plot concerns a part of the DSW (-590 < 
X < —560), and the right one concentrates on the region of 
the obstacle. In both the right and left plots the black solid 
line corresponds to the analytical density profile. In the right 
plot the red solid line is the numerical density computed in 
the case A = 5. The numerical profile is represented by red 
points in the left plot, because it would be indistinguishable 
from the analytical result if represented with a solid line. 



FIG. 6. Density profiles a,t t = 500. The three upper panels 
present the results of numerical simulations performed in the 
case where the obstacle is represented by the Gaussian poten- 
tial ((5T1) and the incident beam by ((50)) with A = 20, 10 and 
5. The lower panel is the analytical result corresponding to a 
point-like obstacle ([9]) and a step- like incident beam (|10[) . (|ll|l . 
i.e., A = 0. 



a train of dark solitons is observed along the plateau in 
region C (cf. Fig. [51 upper panels). As it is clear from 
the figure, the number of solitons in the plateau decreases 
with A and the appearance of a train of solitons is thus 
an effect of the finite width of the front of the initial wave 
packet. These solitons are generated at the initial stage 
of evolution and their precise characteristics depend on 
the details of the initial density distribution. We note 
that the number of solitons in the train is several orders 
of magnitude lower than the number of nonlinear oscil- 
lations in the DSW (see Fig. [6]), i.e., the existence of 
this train of solitons is a minor effect in our asymptotic 
description of the flow. Nonetheless, in any experiment 
the incident wave packet will not be infinitely sharp and 
the resulting train of solitons should be easily observed. 
It is therefore worth spending some time discussing its 
properties. 

The periodic solution ([13]), p4)) describes a train of 
solitons if one considers a situation where A2 — A3 while 
the other Riemann invariants (Ai and A4) remain con- 
stant. Here one has 

^'^'^''JTWc' ^4-Ao. (83) 
The velocity V and minimal density nmin of one of the 



solitons of the train can be expressed as functions of the 
coinciding Riemann invariants X2{x,t) — A3(x,t): 

y=i(Ai+2A3 + A4) = y +A3, (84) 

and 

'^min =x(Al + Xi)"^ - A3(Al + A4 - A3) 

7 (85) 

where ui is given by Eq. (j55|) . Mc being determined as 
a function of k and Aq by Eq. ([54t . In the present case 
(no = 1, "0 = 1 and k = 5.2) one gets ui = 0.039. 
In accordance with the general spirit of our approach, 
we again assume self-similarity of the solution, i.e., the 
locations of the solitons are given by an equation of the 
form X = Vt. The corresponding values of 

A2 = A3 = f-^ (86) 

are shown in Fig.[3]as a (red online) dashed line attached 
to the soliton edge of the DSW Note that dM]) 

corresponds to the solution ([32]) of the Whitham equation 
in the soliton limit where \2{x/t) = \3[x/t). 

Then, elimination of V from ([55]) yields the following 
relation 




between two easily measurable parameters: the minimal 
density rimin of a soliton and its coordinate x at time 
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FIG. 8. The thick solid line represents the analytical predic- 
tion (|87|) for the positions of the points of minimal density of 
the soliton train at time t — 500. The value ui — 0.039 was 
used in the plot of Eq. 1)87^ [see discussion below Eq. (|85|l ]. 
The dots show the same quantity obtained from the density 
profile n{x, t = 500) computed numerically in the case A = 20 
(see Fig. [6l upper panel) . 



t. The very good agreement of this analytical prediction 
with the numerical results is illustrated in Fig. [5] where 
the solid line represents the curve ((87)) whereas the dots 
show the minimal density at the position of the solitons 
as observed in Fig. [6l upper panel (in the case where 
A = 20). 

Hence, although the detailed description of all the 
features of the soliton train — depending of the precise 
shape of the incident front of the wave-packet and of the 
potential — is beyond the scope of our approach, we still 
can get some insight on the loci of the density minima by 
assuming a kind of self similarity which relies on the fact 
— confirmed by inspection of the numerical simulations — 
that the train is formed in a very short period, at the 
initial stage of the collision between the condensate and 
the obstacle. It then evolves according to the solitonic 
limit of Whitham equations. 

VI. CONCLUSIONS 

In the present work we studied the flow of a Bose- 
Einstein beam incident onto an obstacle in a regime 
where a structure forms similar to the supersonic expan- 
sion obtained in a Laval nozzle. For an appropriate choice 



of incident flow and of potential modeling the obstacle we 
showed that, at large time, an acoustic horizon is created 
around the obstacle. Its formation is accompanied by a 
release of energy in the form of a dispersive shock wave 
whose characteristics were studied in detail. This made 
it possible to precisely determine the condition of for- 
mation of the structure studied in the present work (see 
Fig. 13). 

The analytical description has been sustained by nu- 
merical simulations which confirmed our analysis. The 
numerics however revealed an interesting and unexpected 
feature: the ejection of the DSW from the region of the 
horizon is typically accompanied by a train of dark soli- 
tons which are only slowly expelled from the obstacle. 
We argued that this train of solitons is formed in the ini- 
tial stage of collision of the Bose-Einstein beam with the 
obstacle and that its detailed characteristics depend on 
the precise shape of the incident beam and of the obsta- 
cle potential. We also showed that one of its important 
features can be described analytically (namely the line 
of location of the density minima), precisely because this 
train is formed in an early stage. 

The existence of this history-dependent train of soli- 
tons is due to the specific dispersive and nonlinear wave 
mechanics describing the Bose-Einstein condensate. It 
remains to investigate whether this train of solitons does 
or does not hinders the experimental observation of sonic 
analog of Hawking radiation. 
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